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Abstract 

We develop a new three-dimensional multiparticle Monte Carlo (SDmpMC) approach in order 
to study the hopping charge transport in disordered organic molecular media. The approach is 
applied here to study the charge transport across an energetically disordered organic molecular 
heterojunction, known to strongly influence the characteristics of the multilayer devices based on 
thin organic films. The role of energetic disorder and its spatial correlations, known to govern the 
transport in the bulk, are examined here for the bilayer homopolar system where the heterojunc- 
tion represents the bottleneck for the transport. We study the effects of disorder on both sides 
of the heterojunction, the effects of the spatial correlation within each material and among the 
layers. Most importantly, the SDmpMC approach permits us to treat correctly the effects of the 
Coulomb interaction among carriers in the region where the charge accumulation in the device is 
particularly important and the Coulomb interaction most pronounced. The Coulomb interaction 
enhances the current by increasing the electric field at the heterojunction as well as by affecting 
the thermalization of the carriers in front of the barrier. Our MC simulations are supplemented 
by the master equation (ME) calculations in order to build a rather comprehensive picture of the 
hopping transport over the homopolar heterojunction. 

PACS numbers: 73.61.Ph, 73.20.At, 72.20.Ee, 05.10.Ln 
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I. INTRODUCTION 



The physical processes governing the devices based on thin organic amorphous films have 
been the subject of significant interest during the last decade. p| This includes the charge 
transport, exciton creation, and electron-hole recombination in the strong electric field of 
the order of IMV/cm. These and other processes form the basics of functioning of present 
and forthcoming organic electronic devices. In addition, these devices are usually composed 
of more than one organic material. This introduces the organic heterojunctions which often 
play the major role in the device characteristics. This shows both experimentally and in 
some elaborate device model simulations. The structure of the heterojunc- 

tion is known to affect the efficiency, electric characteristics of the device, and its durability. 
Therefore the understanding of the electronic processes at organic heterojunctions is essen- 
tial for the proper understanding of the whole device. In spite of that the theoretical work 
aimed at understanding the organic heterojunctions, in terms of their structure, is scarce. 
This may be contrasted with numerous theoretical investigations on the effect of disorder 
and correlations on injection and charge mobility within the bulk. This former research 

has established a general picture of electronic transport in disordered molecular materials 
and polymers. It proceeds via electron hopping among representative molecular states. In 
particular, one deals with the highest occupied molecular orbital (HOMO) and the lowest un- 
occupied molecular orbital (LUMO). The occupancy of HOMO and LUMO orbitals changes 
as the device is switched on. First, extra electrons and holes are being introduced from elec- 
trodes (injection) onto neighboring molecules. As the carriers spread further into the bulk 
they reach the heterojunctions. The most important characteristic of a hetero junction is 
the difference of LUMO and HOMO levels between different molecular species. This differ- 
ence often shows as the energy barrier for the carriers (electrons entering the material with 
higher LUMO level or holes entering the material with lower HOMO). This almost regularly 
occurs in multilayer structures of organic light emitting diodes (OLED), where lower LUMO 
materials near the cathode (facilitating electron injection) meet the higher LUMO materials 
inside the device (more appropriate for light emission of the photon with several eV's in 
energy). Similarly, one encounters the energy barrier for holes at heterojunctions closer to 
an anode. This stepwise injection, with intermediate materials between the electrode and 
the electro-luminescent material, is known to facilitate the charge inflow into the device. 
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However, it also happens that in that way the bottleneck for charge transport moves from 
the electrodes towards the organic heteroj unctions. The latter then dominate the electric 
characteristics of organic multilayer structures. It should be pointed out from the outset 
that the average value of the energetic barrier at the heteroj unction is not its sole prop- 
erty to determine its permitivity to carrier flow. It is the energetic and structural disorder 
that also affect the performance of the heterojunction. Disorder is expected to modify the 
heterojunction characteristics, similarly as it predominantly determines the carrier mobility 
in an organic amorphous material. Another direct consequence of the barrier is the charge 
accumulation at the heterojunction, which is much more pronounced there than in the bulk. 
Therefore the Coulomb effects are also more pronounced at the heterojunctions than else- 
where. In that respect it is important to emphasize that in organic molecular materials the 
interacting electrons are well localized in space, with the localization length being of the size 
of a single molecule. Thus the nature of the Coulomb interaction is very much different from 
that in band semiconductors, where the spread of the electron wavefunction is considerably 
bigger than the mean carrier separation. The picture of interaction in disordered molecular 
media resembles the one between slowly moving "point charges". Therefore the interaction 
energy depends essentially on the charge configuration at a given instant of time and may 
not be properly expressed through the average occupancy of sites - the effect of correlations 
is pronounced. 

For the reasons given above the most direct theoretical approach to an organic heterojunc- 
tion is to develop a 3D numerical model for the hopping transport of many particles. In that 
respect it may be noted that numerical models for single-particle hopping were successfully 
used in the past to examine the effect of disorder on mobility and injection. 
These models have been conveniently implemented through the Monte Carlo (MC) algo- 
rithm. The MC approach has been used to simulate the hopping movement of a particle 
exploring the energy landscape of the molecular material. In order to consider the hetero- 
junction problem many such particles have to be considered concurrently. 

It should be noted immediately that the rather simple multiparticle MC approach pro- 
posed here is possible because in the systems under consideration the quantum phase of 
the electronic wave function is essentially lost in the (phonon-assisted) hop. Therefore, the 
orbital quantum effects are not expected to show in the assembly of particles. Instead, 
one may consider electrons as classical particles exhibiting hopping transport. On the other 
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hand, the spin state may be, and certainly is, in the absence of spin-orbit interaction centers, 
remembered in the hop. 

There is one additional reason why the multiparticle MC approach is rather appropriate 
for studying organic molecular devices. The approach in which each particle is followed 
through the system profits from the fact that even for a large number of molecules (several 
million in our simulations) the number of particles is relatively low (several hundred). Indeed, 
in real systems like organic molecular diodes the density of particles is usually very low 
because of the lack of countercharges. Therefore the Coulomb forces limit the charge density. 
On the other hand, when both electrons and holes exist within a single organic layer, it is the 
recombination that limits the density. In both cases the molecules greatly outnumber the 
carriers. Only in the regions near bipolar hetero junctions, where both electrons and holes 
accumulate, the density of particles may be significantly increased. However, even then the 
MC approach may turn tractable. 

As it will become apparent, the 3D multiparticle Monte Carlo (SDmpMC) approach 
may be applied to a wide set of problems aside from heterojunctions, e.g. to study the 
recombination crossection, for considering the density and doping effects for the hopping 
transport, or for considering the influence of the discrete charge hopping on the space charge 
effects. As for heterojunctions, it may be applied both to a hetero junction where only 
one type of carriers is present (homopolar) as well as to those where electrons and holes 
meet and recombine (bipolar). Multiplying the effects of disorder, correlations, density 
effects. Coulomb interactions, recombination and interface structures brings us to great 
many number of possibilities, each within the reach of the method. However, we do not 
consider all these issues in the present paper. The paper introduces the SDmpMC method 
and discusses its application to the homopolar hetero junction problem. We consider the 
heterojunctions where only energetic disorder is present, postponing the discussion of the 
structural disorder for a separate paper. With all these limitations imposed it still turns 
out that the problem of a fiat, energetically disordered, homopolar hctcroj unction presents 
a multitude of effects. We believe that all these have to be understood before addressing 
those related to the interface structured in the real space and before considering the bipolar 
problem. 

Confining ourselves to the problem with a single type of carriers we bring in an addi- 
tional simplification - we do not have to take care about spin degree of freedom. In the 
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simulation we do not permit a single molecule to be occupied by more than one carrier of 
a given type (no two extra electrons or two extra holes may reside on the same molecule 
simultaneously), whatever its spin projection. This constraint, imposed in accordance with 
the usual treatment of localized (donor, acceptor) states in doped semiconductors, [l^ is 
maintained within the model even when we switch off the long-range tail (nearest neighbour 
and beyond) of the Coulomb interaction. Effectively, the constraint implies that single state 
is available per molecule for a given carrier. This implies that the spin degree of freedom will 
not show in the calculation for "electrons-only" or "holes-only" devices. It is understood 
that the "spinless-particle" approach may not be appropriate when electrons and holes are 
considered concurrently, with the recombination being dependent on the spin state of the 
particles. 

From here the paper goes as follows. In the next three sections we introduce the Monte 
Carlo model and the numerical procedure as well as the master equation (ME) approach. 
We will discuss first the MC and ME results for the dependence of the particle current on 
disorder strength in variously parameterized heteroj unctions. We then discuss in more detail 
the effect of long-range and short-range Coulomb interaction on the barrier crossing. The 
ratio of the forward and backward hopping rates is discussed in a separate section. This 
ratio also turns out to be important in the final section where we consider in more detail 
the effect of the precise form of the hopping law on the properties of the system. 



II. THE MODEL 



As the model of a heteroj unction we consider the model of two cubic lattices placed side 
by side. For MC simulations we use two lattices of the size 200 x 200 x 50 sites. Each of them 
represents a separate organic layer (denoted as layer I and layer II). The overall thickness of 
the model sample (100 molecular monolayers) is thus typical of real devices. The electrodes 
that are placed on both sides of the device impose the electric field and the carrier injection. 
For the two lateral directions, perpendicular to the applied electric field, we impose periodic 
boundary conditions to account for the extension of the device in those directions. 

Here we describe a somewhat more general model than required for the calculations in 
the rest of the paper. In the present section we consider the system where both electrons and 
holes are present. Each site is then characterized by two energy levels, LUMO and HOMO. 
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The energetic disorder is introduced by adding the random values ej to the site LUMO and 
HOMO energies. These values follow the Gaussian distribution (27ra^)^^/^ exp (— e^/a^). 
We consider both uncorrelated and correlated disorder. For the latter, the short-range 
correlations between the site energies are further imposed by taking the energy of the i-th 
site to be proportional to the average of the energies of all the sites contained within a 
sphere of radius acorr, i-e. Cj = Ncorr^k^k, where Ncorr is a normalization factor, jisl lisj ] 
The correlation in the heteroj unction region is handled separately. We consider both a) the 
case where the correlations inside each device layer are independent, and b) the case where 
the correlations extend over the organic/organic heteroj unction. In the latter case we speak 
of "over-barrier" correlations. This case can occur especially when the energetic disorder 
is caused by random orientation of molecular dipoles. The correlations, occurring due to 
the long range of the dipolar potential, then extend beyond the interface between different 
materials. 

The charge carriers (electrons and holes) interact through the Coulomb interaction, 
qiQj / iTTeQerTij , screened by the dielectric constant of the material. Of course, the "self- 
interaction" of the electron is excluded. In particular this means that the energy levels of 
the two sites involved in hopping have no Coulomb contributions from the electron involved. 
On the other hand, these energy levels do experience the Coulomb shifts from all other 
particles in the device. In practice the calculation of the Coulomb shift is simplified by con- 
sidering only the particles in the basic cell as point charges. For a given site this basic cell 
accounts for the portion of the device, with lateral extension 200 x 200 in our case, centered 

n 

around the particular site. The Monte Carlo calculations for one-component plasma|lfil| 
showed that this approximation (termed "the minimum-image convention"), produces small 
differences with respect to the exact Ewald sum for a diluted plasma. 

Due to weak binding forces in organic molecular semiconductors the carriers are localized 
on molecules. The hopping of electrons and holes from one localized site to another occurs 
when (as a result of the electron-lattice interaction) the atomic vibrations change the relative 
energy of these localized states. Three parameters determine the conditions of the charge 
transfer between the sites: the transfer integral J related to the wave function overlap 
between sites, the coupling strength g of the carrier to the lattice and the typical phonon 



energy Jtwq 
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18l Il9l |20[ In the organic materials used as transport and emission layers in 



OLEDs, the transfer integral J is usually less than 0.1eV.|21,] The phonon mode which 
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produces the relaxation is the C-C stretching mode which interchanges single and double 
bonds of the conjugated backbone {hwo = 0.15 to 0.2eV). The coupling constant can be 
deduced from the relaxation energy of the molecule upon charging: g = EJfvjjQ. Values 
of the relaxation energies have been calculated in Alq3|:22] and in PPV.[22[ They range 
between 0.06 and 0.3 eV. A comprehensive review on the relative importance of the relevant 
parameters and the qualitative effects that we follow here, has been published by Emin. | 
It may be recalled that the ratio J/hjjQ is informative about the adiabaticity of the 
process. In the molecular materials considered here, where J is small, hopping is usually 
nonadiabatic: the intermolecular coupling energy is so small that the electronic carrier can 
only occasionally respond to favorable changes of atomic configurations by moving between 
sites. On the other hand, the coupling strength g in those materials is in the range from 0.5 
to 2, implying a rather strong coupling. The hopping rate for such a process is given bvjiol 



where Ei, is the polaron binding energy. 

However, for historical reasons, hopping in molecular materials is usually considered 
through formula appropriate for a nonadiabatic process in the weak-coupling regime. The 
Miller-Abrahams formula, originally proposed for shallow impurity transport in semicon- 
ductors, gives the jump rate for such a phonon-assisted process. I25I 



where uq is the attempt-to-jump frequency. Here only the rate upward in energy is activated. 

Powered by the success of the Gaussian Disorder Model of Bassler and coworkers 
and due to its simplicity, the Miller- Abrahams formula was used in almost all simulations 

. Nevertheless, in spite of their differences, the two 
formulae lead to some similar physical processes in a disordered system. Field-dependent 
mobility is also qualitatively similar for the two laws. On the other hand, one big difference 
is the appearance in Eq.(0) of the Marcus inverted region. In fact at higher fields, or more 
generally when |ei — ej| ^ E^, the hopping rate starts decreasing due to the decreasing 
probability of phonon emission for a downward hop in energy, jl^ 

In order to provide the link with most of the papers on the hopping transport in the 
systems with Gaussian disorder, in this work we use mainly the hopping processes described 




(1) 
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by Eq.Q. The question of the influence of different hopping laws on our conclusions is 
addressed in the separate section towards the end of the paper. 

III. THE 3D MULTI PARTICLE MONTE CARLO ALGORITHM 

The multiparticle Monte Carlo algorithm we use may be roughly described as follows. 
Basically, we follow the evolution of many particles, electrons and holes. As the system 
evolves through hops, at each step one has to select the particle that is to hop next. For 
that purpose we introduce the "dwelling time" tdwcu- amount of time that the 

particle at site i spends there before the hop, tdweii — ^ij^ ■ "^^^ particle with the 
smallest dwelling time is selected to perform the jump at a given instant of time. Once the 
particle at site k is chosen to perform the next hop, the destination of the hop is determined 
by a random number generator in accordance with the relative magnitudes of the Ukj of sites 
inside a box of size Nb x Nb x Nb (e.g. with Nb = 6). After the hop is performed the 
time variable t is advanced by the corresponding dwelling time - i.e. the time step in the 
simulation is defined by the dwelling time of the particle executing the hop. 

More precisely, the algorithm goes as follows: 

1. Let t denote the absolute time of the simulation and ti the private time of particle 
at site i. The latter is set to zero when the particle arrives at the site and advances 
when the particle do not hop. Initially ti=0 for all sites i populated by particles. The 
total number of occupied sites corresponds to A^, the number of particles (electrons 
and holes) in the system. 

2. Calculate the energy level shifts for the sites inside the hopping box {Nb x Nb x Nb) 
around each particle. Each energy level is determined by the external electric field 
and the Coulomb potential of the particles in the system. 

3. Calculate the hopping probability for each site in the hopping box. 

4. Deduce the dwelling time tj^eii occupied sites. 

5. The site k, occupied by the particle to hop next, is determined by searching for the 
smallest (t^weii ~ ^0 among occupied sites. 
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6. Choose at random the site j within the hopping box of site k as the destination for the 
hop, obeying the relative probabihties calculated above. . The particle is de-assigned 
from site k and assigned to site j. 



7. Update the times: < 



t, = 



ti = ti+ (tdweii - tk) , foroccupiedsites, i ^ j 

If (^dweii ~ tk) < then the particle at site k hops to the site j chosen in the same 
manner, while the time variables are updated as follows: 
t,=0 

t, ti unchanged, for occupied sites i ^ j) 

The new configuration of particles is used to continue the simulation by re-entering the 
above algorithm from step 2. 

The unit of time in the code is given by 

to = [Guoexp {-2'ja)]~^ (3) 

The parameters z/q, and 7 are the same throughout the system. In the simulations presented 
here z/q = 10^^ s~^ and 7 = 8.33 nm^^. 

As the model is currently set for the study of organic/organic interfaces, the details of 
the exact injection mechanism will not matter very much. However, it is important to 
establish the mechanism for the carrier inflow towards the interface. For example, one can 
use an appropriate injection formula with the parameters suitably chosen, e.g. the one 
discussed in Ref. 0|. The probability for the injection at the sites near the electrode is then 
determined by the electric field at the injecting electrodes. The latter is separately calculated 
within the code. Alternatively, one may fix the number of carriers (separately electrons and 
holes) in the device and create a carrier at a random position near the injecting electrode 
whenever a particle disappears either by leaving the device at the opposite electrode or in 
the recombination process. This approach of a fixed number of particles is chosen for the 
present paper, because in that case the total current in the system is the property of the 
heterojunction and not related to electrodes. However, to some extent, the results obtained 
within the two "injection approaches" may be transformed from one into another. For 
heterojunction investigations this basically reduces to calculating the current per particle. 



instead of considering only the value of the current at a fixed number of particles. The 
calculation of the Coulomb effects on the injection field is simplified by the fact that most 
of the carriers are located in the heteroj unction area. 

The simulations related to the present paper begin with a random distribution of carriers. 
Snapshots of all the variables of interest are taken each 50 hops of the fastest carrier (i.e. 
in periods in which at least one of the carriers makes 50 hops). The steady state (the same 
average current along the whole device) is reached within the period of approximately 14 
snapshots. After this stage the average of all the variables is taken in the time window of 
50 snapshots. In this way we obtain satisfactory statistics for typical values of disorder and 
field strengths. 

The algorithm described above has been implemented on a cluster of Linux-based ma- 
chines. The numerical application of the SDmpMC method is much facilitated by the de- 
velopment of the computer cluster calculation techniques, where a handful of particles is 
taken care of by each computer node. The information required to calculate the Coulomb 
interaction between particles is communicated among nodes using Message Passing Interface 
(MPI) calls. 



IV. THE MASTER EQUATION APPROACH 

Many MC methods are known for their conceptual simplicity and a straightforward nu- 
merical approach. However, they are also known for the need to examine great many states 
of a system in order to reduce the statistical error. The present MC method shares this 
property. The other method that is sometimes used to examine the system of particles hop- 
ping in a disordered energy manifold is the master equation (ME) approach. The master 
equation approach has some advantages over the MC approach with respect to reducing the 
statistical uncertainties, as discussed in Ref. [3| 

On the other hand, the ME approach fails to treat the Coulomb interaction correctly, 
both long-range and-short range, although some attempts have been done to consider ap- 
proximately the no-double-occupancy constraint within this approach. 28, Here we use 
this approach where appropriate to supplement our results of the MC simulations. A more 
detailed account on the master equation approach may be found in Refs. 0, 1^ where it 
was applied for studying the mobility and injection in a 3D disordered system, respectively. 
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Here we briefly outline the method and contrast it to the 3DmpMC approach. 

The basic assumption of ME approach is that the occupancy of site i evolves in time 
according to the equation 



where z/j^'s are hopping rates already introduced within the MC approach. For an ensemble 
of non-interacting particles this is a perfectly rightful assumption if rii represents the average 
occupancy of the site over a time scale much bigger than the inverse of the characteristic 
hopping frequency. Finding the steady state of the system {drii/dt — 0) then amounts to 
solving the (huge) linear system of equations for nj's, each of the equations being of the type 
as given above. The source part for the linear systems is provided by suitable additional 
terms that describe the coupling to the particle reservoir (electrodes). The solution simul- 
taneously takes care of all possible hopping processes in the system. This is at odds with 
the MC approach which samples only a fraction of the bonds in the system, and those that 
are sampled are visited a flnite number of times. However, as usual with the Monte Carlo, 
the statistics improves as the sampling extends. 

As already mentioned, the effects of interaction are difficult to introduce correctly within 
the ME approach. For example, although some neighboring sites i and j may be occupied 
to the similar extent on average, nj ~ nj, in reality this may come form the same particle 
visiting both sites in turn. This means that at each moment one site is occupied while 
the other is empty and vice versa. Thus, calculating the interaction energy as the product 
rii X rij is obviously erroneous. 

The simplest interaction that one wants to take account of is the strong repulsion among 
particles on meeting on the single molecule. If this interaction is strong enough the double 
occupancy of a molecule will be practically forbidden. Thus the particle at a site i will not 
be able to jump to a site j if the latter is occupied. In the quantum Hamiltonian language 
this is often described by effectively multiplying the hopping operator by {1 — fij), fij being 
the particle density operator for site j. A similar extension to the master equation approach 
is obviously approximate since Uj in the master equation does not denote an operator but 
an average site occupancy. 




(4) 



11 



However, on the qualitative level the extension of the ME equations to 

seems advantageous. In the first place it assures that, irrespective of the number of carriers 
in the system (which should be less than the number of molecular sites, of course), the 
system will acquire a steady state with < 1. For example, this prevents deep states 
(traps) to be overpopulated when the average number of carriers increases. Moreover, in 
the limiting case of no sources or sinks in the system, the equilibrium will be given by 
rii — {1 + exp[(£^i — n) / kBT]}~^ , with fj, set by the total number of particles. The latter 
result is what one expects for a site with no-double-occupancy constraint being in equilibrium 
with the reservoir represented by the rest of the system. For these reasons the preceding 
generalization of the ME seems suitable for considering a system with a finite concentration 
of particles. On the technical level, the extension is paid off by the fact that the steady-state 
equations become nonlinear. This requires devising a suitable mathematical procedure to 
solve the huge nonlinear system of equations for a given number of particles in the system. 
On the fundamental level, it must be kept in mind that the nonlinear ME approach is 
intrinsically approximate in the sense of the mean-field approximation, as explained above. 
Therefore, the multi-particle Monte Carlo approach, while possibly inferior (statistics-wise) 
for considering a system of non-interacting particles, seems to be an appropriate and superior 
approach for considering a system of real, interacting particles. 

All apphcations of the ME approach in this paper are thus for particles that do not 
interact with long-range Coulomb forces. The linear and nonlinear ME solutions are com- 
pared to assess the importance of short-range repulsion for a given number of particles in 
the system. To reduce the complexity of the computational problem we consider only the 
nearest-neighbor hopping in all ME calculations. The results are generally consistent with 
MC simulations, given the strong decay of i/jj's with increasing distance among molecular 
sites. All ME calculations presented here are related to systems with correlated disorder in 
both layers, the correlation not extending over the heterojunction. 
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V. THE EFFECT OF DISORDER ON THE CURRENT THROUGH THE HET- 
EROJUNCTION 



The first important result of our MC simulations for homopolar lieterojunctions are col- 
lected in Fig. n The figure shows the results for the current in the systems with correlated 
and uncorrelated disorder, with and without long-range Coulomb interactions, and contrasts 
the role of disorder in one vs. both layers in the bilayer device. For layers with correlated 
disorder the correlation length is fixed to three lattice constants. 

We briefly comment on the main features discernible from FigQ deferring a more detailed 
discussion of each of them for separate sections. The figure shows the general dependence 
of the barrier permittivity on disorder strength. The only case where the current increases 
with increasing disorder strength is the one in which disorder is confined to the layer II while 
the layer I is kept ordered. The observed increase of the current with disorder is presumably 
due to the particular points where the barrier across the hetero junction is locally lowered. 
The current over the barrier is generally expected to depend exponentially (i.e. strongly 
nonlinearly) on the barrier height. Hence it is expected that the contribution of those 
crossing points to the average overcomes the effect of the points where the barrier is locally 
increased. 

For all other cases the disorder strength is kept equal in both layers. In all these cases the 
current diminishes as the strength of disorder increases. Thus the effect of disorder in the 
layer I is to reduce the current. This effect always dominates the previously observed effect 
of disorder in the layer II. As discussed later in the paper and already observed in Ref. 
the carriers accumulated in front of the barrier tend to thermalize to low energy states there, 
thus effectively increasing the barrier at the heteroj unction. The figure also suggests that 
correlation within separate layers is not particularly beneficial for barrier crossing. On the 
other hand, the barrier crossing is facilitated when the correlation among layers is present. 
It may be further observed that switching on the long-range Coulomb interaction always 
increases the current in the system. This is chiefly because the Coulomb forces effectively 
lower the barrier. However, as will be shown soon, the size of this effect is much smaller 
than expected from the usual naive argument. 
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VI. THE EFFECT OF THE COULOMB INTERACTION ON THE BARRIER 
CROSSING 



A. Long-range Coulomb forces 

The effect of the interaction may be observed in our simulation as we increase the number 
of particles in the system. This is illustrated in Fig|21 The figure shows the current as the 
function of the number of particles in the system for particular values of the barrier height 
and disorder strength. The points related to the system of particles without long-range 
Coulomb interaction follows a straight line as the number of particles increases. The effect 
of the long-range Coulomb interaction is unobservable for small number of particles in the 
system. As expected, the effect starts to show as the number of particles increases. 

It is also expected that long-range Coulomb interaction among electrons accumulated 
at the hetero junction will help the carriers to pass the barrier. This is expected since the 
electric field due to the interaction may lower the barrier in a way similar to how the external 
electric field F lowers the barrier by the amount Fqa. The usual straightforward estimate 
for the barrier lowering due to particle interaction goes as follows. Let us denote the average 
planar charge density at the hetero junction by ps = qN^/a'^, where A^^ stands for the average 
number of carriers per cell cross section a^. If the charge at the heteroj unction is assumed to 
be smeared homogeneously over the heteroj unction plane, the extra electric field produced 
by this charge would be Fhom = Psf^i^o^r)- Twice that much is the change of the electric 
field over the heterojunction. The naive estimate for barrier lowering is Fi^omQo,- 

However, it must be recalled that the carrier wavefunction is not at all homogeneously 
smeared in the heterojunction plane. The carriers localized on molecular sites are spaced by 
d ~ ^Jl/Ng on average. It is the Coulomb repulsion among accumulated electrons that tends 
to order them in a regular mesh in the heterojunction plane. The continuum approximation 
-^hom for the electric field is valid only for distances 5x ^ d from the heterojunction. At 
shorter distances the effect of charge discreteness must be accounted for. Specifically, the 
distance to be considered for barrier lowering is 5x ~ a. As the characteristic case for a 
homopolar heterojunction is a <^ (i, the total change of the electric field at the heterojunction 
is expected to be smaller than the average field in the smeared charge picture. The effective 
electric field responsible for barrier lowering may be estimated e.g. by assuming that repelling 
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charges order approximately in the regular square lattice. The effective electric field then 
turns out to be very much smaller than Fhom, as showed in Fig. El 

We illustrate these considerations with the case of 140 particles in our 100 x 200 x 200 
system. The particles accumulate in front of the heterojunction, as shown in Fig. 0] The 
planar particle density per unit cell cross section is 140/(200 x 200) ~ 0.002, with most 
of the particles concentrated in the last few monolayers of the layer I. The surface charge 
density is ps = 140g/(200a)^. In the continuum approximation the accumulated charge 
produces the extra field of approximately 0.25 MV/cm (taking = 3.5). This field is 
already quite strong. Twice as many particles would double the electric field in layer II 
with respect to the imposed 0.5 MV/cm and reduce the electric field in the layer I to 
zero. For that reason the number of particles in the system may hardly go much above this 
concentration. Further particles would cause a negligibly small or even negative electric field 
at the injecting electrode. Taking the calculated extra field of 0.25 MV/cm as responsible 
for barrier lowering would imply a barrier shift of 0.015 eV. This would imply a current 
amplification of approximately (exp(0.015eV/A;BT) ^ 1.8) 80 percent. However this is not 
what one observes in the MC simulation where the amplification of a few percent is observed. 

On the other hand, taking into account the discrete nature of the charge distribution in 
the heterojunction plane gives a much smaller electric field. It turns to be approximately 
20 times smaller than the one obtained for the homogeneously charged plane. The barrier 
reduction due to particle interaction is much smaller than obtained above, as well as the 
current amplification which is then estimated to (exp [(0.015eV/20)/A;B7'] ~ 1.03) a few 
percent, as observed in the simulation. 

In summary, the Coulomb effect usually calculated for a homopolar heterojunction is 
strongly overestimated. This is true for most multilayer device models, as well as for some 
recent analyses of the dominant device physics. 0] On the other hand, the effect of discrete- 
ness that reduces the Coulomb effects in homopolar heteroj unctions is expected to strongly 
amplify the Coulomb effects among electrons and holes at bipolar junctions. 

B. The hard-ball repulsion 

The short-range repulsion is always present in our SDmpMC model. Therefore its effect 
can not be deduced directly from the results on the current density. On the other hand, we 
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did the calculation with and without no-double-occupancy term within the ME approach. 
We found that the current enhancement due to short-range repulsion at a given particle 
density is less than one percent for the parameters given in FigC] 

The effects of the short-range repulsion can be assessed even within the MC approach 
by considering the distribution in energy of the particles accumulated in front of the energy 
barrier. The particles in the sample spend much more time in front of the barrier then 
elsewhere. The number of hops per unit time among sites within that monolayer greatly 
outnumbers those that go down-field. This makes it possible for particles to explore the 
sites in front of the barrier rather well and establish a local quasi-equilibrium. 

The properties of the quasi-equilibrium with and without the hard-ball repulsion may 
be approximately calculated in a simplified manner, and the results may be compared 
against the MC results. The calculation goes as follows. For non- interacting (n^) par- 
ticles, the probability of a site being occupied is given by the Boltzmann expression 
fsiE) = exp{—E/kBT + ^ni/kBT), where fini stands for the local chemical potential, which 
depends on local concentration of particles. The distribution of the sites in energy, around 
zero average energy, being given by Po{E) = (27ra)-i/2gxp(_E72(T2), the probability dis- 
tribution for the energy E representing an occupied site goes as 



= PoiE)fB{E) = -==exp 



V2TTa 



2 



+ 



2 ^ 2^2 keT 2A;|T2 



(6) 



Thus the mean energy of the occupied state is shifted by —a^/ksT with respect to the 
case without disorder. This well-known shift of the energy of occupied states due to 
thermalizationjll| tends to increase the effective barrier at the heteroj unction. 

This shift decreases to some extent once the hard-ball [hb) repulsion is taken into account. 
The no-double-occupancy constraint changes the probability of the state of energy E being 
occupied to fhb{E) = {1 + exp[{E — fihb)/kBT]}~^ , i.e. to the Fermi-Dirac distribution 
function. ^32] • The chemical potential and mean energy have to be calculated numerically 
for this probability distribution. The result for two cases (i.e. with and without hard-ball 
repulsion) is illustrated in Fig. In general, the hard-ball repulsion reduces the shift of 
the mean particle energy due to disorder. The interaction also shifts the chemical potential 
upwards. The MC results are also shown in Fig. El The thermalisation obtained within 
MC simulations is somewhat weaker than the one predicted by simple hard-ball statistics 
calculation. We believe that this discrepancy may be due to the confining effect of the electric 
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field. By pushing the carriers towards the interface, the field prevents them from exploring 
a larger energy manifold with lower energy levels. By switching on and off the long-range 
Coulomb interaction we checked that those contribute negligibly to thermalisation at a given 
particle density. 

Fig. also shows that for a particular particle density and for room temperature, the 
difference between the thermalization shifts is very small for disorder strengths less than 
0.06 eV. The difference becomes pronounced in systems with stronger disorder or at higher 
particle density at the heterojunction. Since the value for a in some disordered organic 
molecular material may be rather big (e.g. a ~ 0.15 eV for Alqs), [s^, the effect of short- 
range repulsion is realistic to expect for some heteroj unctions. On the other hand, the 
particle density at homopolar heterojunctions may be much higher only if the operational 
field is raised much above 0.5 MV/cm. This is because, the operational field limits the 
particle surface density at the heterojunction, as already pointed above. 



VII. BACKWARD AND FORWARD HOPPING RATES 

Recently Arkhipov et al. got interested in the ratio of backward (6) to forward (/) 
hopping rates at the heteroj unction [3]. Their physical motivation was to resume the hops 
that lead to injection at a heterojunction (oc 6 — /) from those that may lead to immediate 
recombination with a carrier of the opposite polarity. While their argument may certainly 
need to be amended by considering the Coulomb interaction among carriers of opposite 
polarity in the latter case, the h/ f ratio is indeed important when considering the permitivity 
of the barrier to current flow. 

Other motivation to consider the h/ f ratio comes from reflecting on the aging processes 
in organic devices. A decade ago it was shown by Tang and Van Slyke that the aging rate 
of OLED's is proportional to the operating current. j^l This rule, also termed the "Coulomb 
aging" rule, may suggest that the hopping rate at the microscopic scale is directly related 
to degradation [3J|. However, while the average current density is constant throughout the 
device, it is not necessarily so for the hopping rate. This is because at any given point in the 
device both forward and backward hops contribute to the local hopping rate (oc f + b), while 
the current, being constant along the device, is proportional to their difference (oc f — b). 
The hopping activity, here defined as {f + b) / {f — b) = (1 + 6//)/(l — b/ f) is therefore 
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expected to be dependent on position in an inhomogeneous device, and especially when 
considering the region near the heterojunction. FiglHl shows this dependence as obtained 
in our ME calculation. It may be observed that the hopping activity is many times bigger 
near the heterojunction than elsewhere. However, it is not the very heterojunction where 
the hopping activity is maximal. The hopping activity maximizes in the last few monolayers 
of the layer I, just before the barrier. This is mostly due to the increased concentration of 
particles in that region and the relatively low probability for further down-field hops. 

Concentrating on the hops over the heterojunction, which determine the heterojunction 
permittivity to particle flow, we investigated the effect of the structure of the heterojunction 
and the Coulomb interaction on the h/ f ratio. The result is shown in Fig. [3 While generally 
the b/ f ratio increases with the increasing strength of disorder, one may also notice the effect 
of correlations and the Coulomb interaction. First, correlation generally lowers the b/ f ratio, 
especially as disorder increases. The effect is possibly related to the fact that the carriers that 
enter the layer II face a locally more ordered environment, which implies smaller probability 
for returning hops. Second, the Coulomb interaction increases the relative importance of 
backward jumps to some extent, especially at low disorder strength. We do not have a 
simple explanation for the latter effect. Considering the forward hopping rate separately in 
Fig. ISl we find that it increases on switching on the Coulomb interaction, to the effect of 
increasing the current, in spite of an increased b/f ratio, J (x f x (1 — 6//). However, the 
increase of the b/ f ratio cannot be explained simply in terms of increased effective electric 
field at the interface. This will become clear shortly as we consider the dependence of b/f 
ratio on the electric field. 

For the nearest-neighbor hopping and without long-range Coulomb interaction, an ana- 
lytic expression for the b/f ratio can be derived in the limit of zero disorder. The derivation 
starts from the fact that the density in an ordered system is constant throughout layer II. 
This establishes the relation between the current and the particle density in the same layer. 
The current may be also expressed by the difference between forward and backward hopping 
rates at the very heterojunction. Then, assuming A > Fqa, the result for the M.-A. hopping 
formula follows, 

(6//)hj = 1/ [2 - expi-Fqa/knT)] (7) 

The zero-disorder value of b/f is 0.59 for the parameters of Fig. [71 in accordance with MC 
results. 
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As for the dependence on disorder strength cr, the MC results are supplemented here 
by the ME results for h/ f away from the heterojunction (FigE)). The dependence on the 
electric field given in FiglTIH while roughly following the analytical result for cr = 0, changes 
differently with disorder for low and high electric field, the crossover being roughly defined 
by Fqa ~ A/3. 



VIII. THE EFFECT OF THE HOPPING LAW ON THE BARRIER CROSSING 



One of the questions that arise when discussing the effects in a system with hopping 
conduction is the qualitative dependence of these effects on the precise form of the hopping 
law. While the results presented above were obtained by assuming the Miller-Abrahams 
hopping law, qualitatively similar results are obtained if one uses the symmetric hopping 
law or the small-polaron hopping formula. This is illustrated in Fig^2 (ME calculations). 
The figure shows the current through the heterojunction as a function of the barrier height 
A, at a given disorder strength. The hopping laws used for this calculation are 

Ei — Ej + \Ei — Ej\ 



.MA 



SYM 



POL 



Unn exp 
J^nn exp 
Unn exp 



Ej — Ej 
2kBT 



Ei — Ej {Ei — Ej] 



(8) 
(9) 
(10) 



2kBT SkBTEb 

To simplify the comparison, the bare value for the nearest-neighbor hopping frequency z/„„ 
was chosen the same for all three hopping laws, while the hopping was confined to the nearest 
neighbors. The number of carriers in the system was also kept the same. 

For A = the current through the system reflects the value of the mobility in a ho- 
mogeneous system for different hopping laws. The issue of mobility vs. hopping law in a 
disordered system has been addressed in the literature. 35| It was shown numerically that 
the Miller-Abrahams and the symmetric formula give qualitatively the same field depen- 
dence of mobility, as far as the disorder strength is bigger than the potential drop among 
neighboring molecules due to the electric field, cr > Fqa. While sharing the common field 
dependence, the value for mobility ^ma{F) in that region is found generally smaller than 
f^SYM{E). We are not aware of any similar considerations for heterojunctions, either numer- 
ical or analytical. Fig. ITT] shows the dependence of the current through the junction on the 
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barrier height. It shows that the dependence is similar for all three hopping laws. Although 
the polaron binding term reduces the current somewhat as the barrier height increases, the 
dependences shown by all three curves are qualitatively similar. 

What may seem particularly strange is the fact that the ratio of Jsym/Jma changes 
very little over the range where the factor exp {A/2kBT), relating different hopping laws at 
the heteroj unction, changes over three and a half orders of magnitude. This may be even 
stranger when the distributions of the carriers in space are examined for the three hopping 
laws. No significant differences can be observed - almost all the carriers are grouped in the 
first few monolayers in front of the barrier, as illustrated in FiglU Therefore the forward 
hopping rate at the heterojunction is approximately exp [(A — Fqa) /2kBT] times bigger for 
the symmetric hopping law then for the Miller-Abrahams (MA) hopping law. Hence the 
backward hopping rate in the symmetric case has to compensate the forward hoping rate 
rather precisely in order to bring the current to the same order of magnitude as for the MA 
case. 

Fortunately, these considerations may be supported by exact analytic expressions for a 
limiting case of the heterojunction between two ordered layers. As already pointed out in 
the previous section, in that case the ratio of backward to forward hopping rates may be 
derived analytically. jSGj For the Miller- Abrahams hopping rate the result was given already 
in EqI7| whereas for the symmetric law it is 

HJ,SYM ^ 1 - (e^W/c.T _ 1) ' (11) 

Neglecting the minor differences between carrier distributions, the ratio of currents in the 
two cases becomes 

JsYM ^ (A-Fga)/2fcBT l ~ iP/ f)HJ,SYM ^^2) 

Jma 1 - (&//) 

The analytical expressions for 6//'s can be substituted into this equation. It is then readily 

verified that the big difference in forward hop frequencies between the two cases is com- 
pensated by respective backward hop rates. This analytic expression for Jsym / Jma also 
reproduces approximately the value obtained numerically (FiglTT|l for finite disorder. We 
believe that this qualitative result also stands in the core of the aforementioned similarity 
between mobilities for the two hopping laws in homogeneously disordered systems. In the 
latter case the bottlenecks for transport are again barriers bigger then Fqa, statistically gen- 
erated and scattered in space. The argument relating the two hopping laws is still expected 
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to apply on those barriers. 



IX. CONCLUSION 

In summary, we presented a new multiparticle Monte Carlo method that is suitable for 
addressing a number of presently unexplored questions related to electronic transport in 
organic disordered materials. In this publication we applied the method to a homopolar 
hetero junction, considering the influence of disorder and Coulomb interaction on the cur- 
rent through the hetero junction. Supplementing our MC results with the master-equation 
(ME) calculations, we developed a rather detailed insight into the effects of short-range and 
long-range Coulomb repulsion, forward vs. backward hopping processes and the precise form 
of the hopping law. In particular, we showed the estimates based on a continuum approxima- 
tion strongly overestimate the effect of the long-range Coulomb interactions on the current 
through the hetero junction. Also, the long-range Coulomb interactions do not affect visibly 
the thermalisation of carriers in front of the barrier. It is rather the no-double-occupancy 
constraint which dominant ly determines the effective barrier height. The presence of disorder 
and especially the energy barrier at the hetero junction lead to position-dependent forward 
and backward hopping rates. The frequency of charging and discharging of molecules then 
varies throughout the device. This likely affects the aging and the degradation. As in the 
case of bulk transport, the microscopic hopping law appeared to play a minor role in het- 
ero junction barrier crossing provided that the polaron binding energy is smaller than the 
energy barrier. This study presents basis for further theoretical studies of more complicated 
(mixed, graded, rough) homopolar heteroj unctions that are being developed experimentally. 
The multiparticle MC method introduced in this paper is currently applied for studying the 
recombination processes at bipolar heteroj unctions. 
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Figure captions 

Figure 1: (MC) Current across the hetero junction as a function of disorder strength. 
Following cases are considered: no corl+ Coulomb: non correlated disorder and (long- 
range) Coulomb interactions included; corl+ Coulomb: correlated disorder within each 
layer and Coulomb interactions included; corl- Coulomb: correlated disorder within 
each layer and without long-range Coulomb interaction; one-side corl-h Coulomb: cor- 
related disorder in the second layer while the first layer is ordered, the Coulomb 
interactions are included; over-barrier correlated -hCoulomb: the correlations in the 
energy disorder present in both layers extending over the hetero junction, long-range 
Coulomb forces being included in the calculation. The number of particles in the 
system (100 x 200 x 200 sites) is always 140. The bare barrier height (the difference 
among average LUMO/HOMO energy of two materials) at the hetero junction is 0.15 
eV. The strength of externally apphed electric field is 0.5 MV/cm. 

Figure 2: (MC) The dependence of the current in the system on the number of 
particles in the system (100 x 200 x 200 sites). The squares stands for the results 
where the interaction among particles is full Coulomb interaction. The lower curve 
(filled circles) stands the current in the system where only "hard-ball" short-range 
interaction is taken into account (i.e. double occupancy of a site is forbidden but no 
long-range repulsion is accounted for). The imposed electric field is 0.5 MV/cm. The 
bare barrier height is A = 0.15 eV. 

Figure 3: The effective field at the homopolar hetero junctions leading to barrier 

lowering, as produced by the accumulated charges, vs. the electric field of the equiv- 
alent homogenously smeared planar charge. Fgfj is estimated here by assuming that 
electrons at the hetero junction are ordered as point charges forming a regular square 
mesh. 

Figure 4: (MC) The average occupancy of the molecular site in the system as a 
function of the position of the site with respect to the plane of the hetero junction (the 
first monolayer of the layer II is positioned at zero) . The density in front of the barrier 
is much bigger than in the bulk of both layers. Prom the heterojunction to the left the 
occupancy falls approximately exponentially. The model parameters are: a — 0.04eV, 
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A = O.lSeV, a — 0.6nm, F — 0.5MV/cm. The data represent the results for for 140 
particles in the system of 100 x 200 x 200 sites. 

Figure 5: The energy shift due to thermahsation in front of the barrier, shown as 
the function of disorder strength. Boltz. shift is for non-interacting particles when the 
Boltzman law gives the orbital occupancy. HB shift is the calculation for the particles 
with hard-ball interaction in the equilibrium for the average molecular occupancy of 
2 X 10~^. The filled circles is what we get from the MC simulations. 

Figure 6: (ME) The hopping activity in the system is shown as the function of 
the position along the system. Here / and b stand for the forward and backward 
hopping frequency between successive monolayers, respectively, f + b is the average 
frequency of hops in both direction, while / — 6 is constant throughout the system, 
being proportional to the particle current in the field direction. Three graphs corre- 
spond to three hopping rules (MA: Miller- Abrahams, SYM: symmetric hopping law, 
POL: small polaron hopping formula with polaron binding energy = 0.05 eV). 
a = 0.05 eV, a = 0.6 nm, A = 0.15 eV. Backward hops at the heteroj unction are 
the most pronounced for the symmetric hopping law, while less pronounced for the 
Miller- Abrahams and small-polaron hopping formulae. 

Figure 7: (MC) Variation of the ratio of backward to forward hopping rates with 
disorder. F — 0.5 V/cm. A — 0.15 eV. The number of particles in the system is 140. 

Figure 8: (MC) Number of forward hops as a function of disorder. The forward 
hopping rate is strongly increased by long-range Coulomb forces, especially at low 
disorder. It is also increased by the disorder correlations between layers ( "over barrier 
corl"). 

Figure 9: (ME) The ratio of the backward to forward hop frequency as the function 
of the strength of disorder. Different graphs refer to different positions with respect 
to the barrier, as indicated by the labels (e.g. -3 denotes the position that is three 
monolayers before the barrier). The spatial dependence of 6// is not at all pronounced 
in layer II. Throughout the layer II the b{a)/ f{a) curve resembles the one labeled by 
-8, corresponding to the monolayer of layer I that is rather far from the heterojunction. 
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Figure 10: (ME) The ratio of backward and forward hopping frequency at the very 
barrier as a function of the electric field at the barrier. The dashed line represents the 
result without disorder, described by the formula in the text. 

Figure 11: (ME) The dependence of the current through the junction as the function 
of the barrier height. The three curves represent the results for three different hopping 
laws. The density of particles per device cross section, 140/ (200*0. 6nm)^ ^ O.Olnm"^, 
is the same in all three cases. 
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